Low-temperature density matrix renormalization group using regulated polynomial 

expansion 
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We propose a density matrix renormalization group (DMRG) technique at finite temperatures. 
As is the case of the ground state DMRG, we use a single-target state that is calculated by mak- 
ing use of a regulated polynomial expansion. Both static and dynamical quantities are obtained 
after a random-sampling and averaging procedure. We apply this technique to the one-dimensional 
Hubbard model at half filling and find that this gives excellent results at low temperatures. 

PACS numbers: 02.60.Cb, 71.10.Fd, 05.70.-a 



The density matrix renormalization group (DMRG) 
method^ is a powerful numerical technique to investigate 
various properties of low-dimensional strongly correlated 
electron systems. The ground state properties are accu- 
rately calculated through targeting the ground state in 
each step of the DMRG process. For dynamical quan- 
tities, a multitarget procedure has been proposed 2 - and 
provides accurate description of various excitation spec- 
tra. These successes of DMRG come from proper choices 
of the target state to be necessary to construct the den- 
sity matrix that contains important bases for the precise 
description of physical quantities. 

The extension of such a targeting procedure to ther- 
modynamic properties has been done for one-dimensional 
(ID) spin systems^ Several tens of lowest-energy eigen- 
states of the systems have been taken as the target 
states and the density matrix is constructed by weight- 
ing a Boltzmann factor for each eigenstate. Later, the 
transfer-matrix method^ has been introduced to calcu- 
late thermodynamic properties in the infinite-size sys- 
tem. The DMRG method is employed for the calcula- 
tion of the maximum eigenvalue of a transfer matrix that 
gives the free energy of the system. Recently another 
finite-temperature DMRG method has been proposed as 
a generalization of time-dependent DMRG£ These finite- 
temperature DMRG techniques successfully give excel- 
lent results for ID spin and electronic systems. However, 
there are some difficulties in each technique; for exam- 
ple, the transfer matrix method is not easily applied to 
complicated models since the transfer matrix is not Her- 
mitian. In the time-evolution method, low-temperature 
properties are not easily obtained since long time evolu- 
tion of the system is necessary. Therefore, it is desired 
to develop a variety of finite-temperature DMRG tech- 
niques, from which one can choose the best one suitable 
for a given model. 

In this Brief Report, we propose a scheme of DMRG at 
finite temperatures, which is a straightforward extension 
of the target-state procedure at zero temperature. The 
target state is weighted by a Boltzmann factor. Making 
use of the polynomial expansion and random sampling, 
we can calculate static and dynamical quantities at finite 
temperatures. In order to obtain good convergency at 



high temperature, we need a large truncation number of 
the density matrix. The proposed method is, therefore, 
suitable for lower temperature region. As a demonstra- 
tion of the method, we show the specific heat, spin-spin 
correlation function, and dynamical current-current cor- 
relation function of the ID Hubbard model at half filling. 
The DMRG results reproduce the exact digitalization re- 
sults at low temperature. 

The DMRG procedure at zero temperature requires a 
target state in order to obtain the ground-state proper- 
ties. Even for finite temperatures, it may be possible to 
have a target state suitable for the evaluation of physical 
quantities. A possible target state may be given by 
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where H is the Hamiltonian, |£) is a normalized arbi- 
trary vector, P is the inverse temperature 1/T, N is the 
dimension of the superblock, and a n = (e n |£), with |e„) 
being the eigenvector corresponding to the eigenvalue e n . 
The inner product of Eq. ([T|) gives the partition function 

Z, provided that a 2 = 1: Z = (£|f) = En=i e ~ 0£n - 
Therefore, Eq. |T]) is a good candidate for the target in 
the DMRG procedure. However, it is difficult to obtain 
all of the eigenstates |e n ) for the superblock Hamiltonian 
whose size of the Hilbert space is of the order of 16m 2 
in the case of the single-band Hubbard model, with m 
being the truncation number of the density matrix. We 
thus need to develop a different technique to treat the 
operator erP H l 2 precisely without obtaining |e„). 

By using the Legendre polynomial expansion for the 
delta function, i.e., 
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with wi = 2/(21 + 1), the Boltzmann factor reads 



-0E n 



(2) 



(3) 



1=0 



where E n is an eigenvalue rescaled to be confined within 
the interval of [— 1, 1]. The corresponding rescaled Hamil- 
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tonian H s is defined as H s — wh(H — A) with scaling 
parameters wh and A. 

In general, there appear so-called Gibbs oscillations in 
any polynomial expansion including the Chebyshev poly- 
nomial often used in the literatures^ The oscillations 
can be eliminated by introducing the Gaussian distri- 
bution function for Pi(E n ) in Eq. $Z§& The polynomial 
regulated by the Gaussian is defined as 
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where a is the half width of the Gaussian distribution 
function set to be 2tt/L, where L denotes the highest 
number of I in the expansion. Inserting the Boltzmann 
factor [Eq. |[3])] into the target state [Eq. (Q])] and return- 
ing to the operator representation, we obtain 

10 * r dee-^ 2 J2 w i lp ^( P ^s)) t r\0- (5) 

Since the integration in Eq. §5§ with respect to e leads 
to the modified spherical Bessel function ii (x) of the first 
kind, the target state is finally written as 

L 

||) ~ C(/?)^ r V-/V2)<fl(£ s )>a 10 , (6) 
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where C(/3) is a normalization constant. 

In order to calculate (Pi(H s )) a we employ a coali- 
tional recursive relation^ 

2/ 4- 1 - 

(Pi+i(H s )) a \0 = J ± T H a {Pi(H a )) tr If) 
-^(P-i^)), 10 + ^a 2 {P((H s )) a |0 ,(7) 
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+(PU(Hs))*\0, 
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where P/(e) ee dP t (e)/de. Starting from (P (H s )) a |0 = 
|0 and (Pi(H s )) a |0 — 10; we recursively calculate 
{Pi{H s )) a |0 U P to / = L and construct the target state 
in Eq. ©. 

In the DMRG procedure, physical quantities are mea- 
sured when the system size is reached to a given num- 
ber in the infinite-size algorithm or enough convergency 
is obtained in the finite-size algorithm^ At this stage, 
we need to introduce a technique to guarantee the re- 
lation o^ = l for the coefficients in Eq. (J7XJ) - This is 
achieved by taking the random sampling of the state |0 
and averaging over the samplings. Let us represent a 
randomly generated |0 as |0 = J2i r i\£i)> where |&) 
is the basis state of the system and is a normalized 
random number generated from a rectangular distribu- 
tion whose center is at zero. Expanding the eigenstate 



|e n ) also in terms of [&), i.e., |e n ) = J2i K,i we ob- 
tain a 2 = Y,i r i b l,i + 2Y,i 7 tj r i r j b n,ib n ,j- After averag- 
ing over many samplings whose number is M, rf will 
become a constant approximately independent of i, and 
riTj will vanish according to l/^/m 2 M^ Therefore, a 
relation a 2 a 1 (n = 1, • • • , N) is expected to be satisfied. 

Physical quantities that do not commute with H are 
also obtained by using the same random sampling. An 
expectation value of an operator A is given by 

(f|i|f) = £4 e -^»<e n |i|e„> 

n 

+ J2 a n ame-^ +em)/2 (e m \A\e n ) . (9) 

The coefficient for the off-diagonal term is expressed as 
a n a m = J2i r i b n,ib m ,i + ;nrjbn,ibm,j- Since rf is 
a constant and r%rj is zero after the sample averaging, 
a n a m is expected to be zero. This means that Eq. © 
gives thermodynamical average of a given quantity. 

We can also calculate dynamical quantities at finite- 
temperature. A dynamical correlation function for an 
operator A may be defined as 
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where 7 is a small positive number. In order to obtain 
Xa (w), we introduce a following expression: 
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with |e) = Y% =0 Wi 1 Pi(e)(Pt(H,))<, \$, and e s = e/w H + 
A. If we take L — > 00 and use the random averaging, we 
can easily find that xa(v) is identical to xa(w)- We use 
Eq. (fTTj) at the stage of the measurement. 

Before reaching the measurement, we need to perform 
a dynamical DMRG procedure using a multi-targeting 
technique i 2 -^ In order to make use of |0 as one of the 
multi targets, we introduce alternative expression 
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which gives xa(^) after taking L — > 00 and the random 
averaging. In addition to possible target states may 
be A |0 and j\ de[w -H + e s - i^ie"^/ 2 |e). How- 
ever, the later state is not easily calculated since it con- 
tains the integration in terms of e. Instead of this state, 
we introduce a state with simple form [u> — H + E — 
i7] _1 A|0, by replacing e s to an e-independent quantity 
E = (£|7f|0- This replacement is based on the fact that 
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the dominant contribution of e s comes from the energy 
range where the product of the density of the eingen- 
states and the Boltzmann distribution function is large. 
In spite of such a rough approximation, we will find that 
this works practically well as discussed below. If this 
replacement does not work well, we should employ the 
original state with the integration as a target state. 

In the present DMRG approach, there are three pa- 
rameters: the number of sampling M, the polynomial 
expansion truncation number L, and the DMRG trunca- 
tion number m. Among them, M is dependent on phys- 
ical quantities as will be mentioned below. The number 
of L in Eq. © predominantly depends on temperature 
T. The lower T is, the larger L is. However, it is diffi- 
cult to obtain numerically the several hundred order of 
il(—/3/2). In such a low-temperature (large (3) region 
where large L is required, we introduce a smaller /3' with 
a relation /3 = nj3' (n is a positive integer), and then per- 
form the polynomial expansion in Eq. ^ n times starting 
from |£) at f3' and inputting the obtained |£) into At 
low temperature, although large L is required, m can be 
reduced as compared with that for higher temperature, 
since the number of the basis necessary to describe low- 
temperature properties is small. As a result, the com- 
puting time is shorter at low temperature than at high 
temperature, in order to get the same level of conver- 
gence. 

We apply a different finite-temperature DMRG 
method to the ID Hubbard model at half filling to 
check its efficiency. The Hamiltonian is given by H — 
~t Ei )(7 ( c U Ci + 1 - ff +h - c -) + U E t ni,tn id , where c[ a (c i>CT ) 
is a creation (annihilation) operator of an electron with 
site i and spin a, n^ a — cj a Ci^, t is the hopping integral, 
and U is the on-site Coulomb repulsion. We use a lattice 
with open boundary condition. 

Figure 1 shows the specific heat obtained by using the 

formula C V (T) = (N S T 2 )- X ((£\H 2 \£) - where 
N s is the number of sites. The result for N s — 8 and 
U/t = 8 is shown in Fig. 1(a), where the exact C v is 
also plotted for comparison. We employ parameters of 
m = 50 and 64 , L = 80 (T/t = 0.1), and M = 1000. 
The error bar due to the sampling is within the size of 
the symbols. At m = 50, the low-temperature peak of C v 
agrees with the exact result. However, with increasing T 
the deviation from the exact result is enhanced. This is 
improved if we use larger m. In fact, taking m = 64, we 
get complete agreement between the exact and DMRG 
results, since 16m 2 exceeds the dimension of the Hilbert 
space. Such an improvement is also clear in the case of 
N s = 20 as shown in Fig. 1(b). Temperatures are re- 
stricted to a range just below the first peak of C v . The 
exact result in the thermodynamic limit iV s — ► oo (Ref. 
12)is also shown for comparison. We find that m = 200 
is enough to get good convergence below T/t = 0.12. At 
higher temperatures, a larger m is required but it is still 
accessible by using a standard computing system. From 
these results, it is apparent that this DMRG method is 
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FIG. 1: Specific heat as a function of temperature in the ID 
Hubbard model at half filling with U/t = 8. (a) iV s = 8 and 
M — 1000. Dotted line shows the exact result obtained by 
the direct diagonalization. (b) iV s = 20 and M = 400. Dash 
line shows the exact result in thermodynamic limit. (Ref. 12) 



efficient for the calculation of low-temperature proper- 
ties. 

As a quantity whose operator does not commute with 
if, we choose the spin-spin correlation function S(r) = 
(£\Sf Sf +r \£) , where Sf is the z component of the total 
spin operator at site i. We choose the two sites, i and 
i + r, in order to make the central site of a given lattice 
the middle of them. Figure 2 shows the staggered cor- 
relation S(r)(—l) r at several temperatures for iV s = 20 
and U/t = 8. The truncation number m is changed with 
temperature in order to get good convergence. The spin 
correlation decreases with increasing temperature as ex- 
pected. 

Finally we show the dynamical current-current corre- 
lation function in Fig. 3. The operator A in Eq. (fTTj) is 
replaced by the current operator j — itj^i <r( c l a c i+\,cr ~ 
h.c). In this calculation, we employ 7 = 0.2i. Although 
7 — > is desired in general, the finite value is introduced 
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FIG. 2: Spin-spin correlation function in the 20-site Hubbard 
chain at half filling with U/t = 8. M = 400. 
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FIG. 3: Dynamical current-current correlation functions in 
the eight-site Hubbard chain at half filling with U/t = 10. 
M = 16. 



here in order to reduce the computational time, in par- 
ticular, at high energy region of Xj(uj). We compare the 
DMRG results of Xj with the exact one for N s — 8. At 
T/t = 0.1, we obtain the same result as the exact one 
even for m = 50. At a higher temperature T/t = 2, how- 
ever, agreement is less satisfactory, since the number of 
m is not enough. This result again demonstrates that 
this DMRG technique works well, in particular, at low 
temperatures. We note that the position of the lowest- 
energy peak dose not change with increasing T, which is 
a consequence of the spin-charge separation^ 



We have shown that it is necessary to perform random 
samplings of the state |£) to calculate physical quanti- 
ties in the process of the measurement in DMRG. The 
number of M necessary to get small statistical error is 
denoted in the captions of each figure. We find that M 
is dependent on physical quantities. We, thus, need to 
check an adequate M for each quantity by examining the 
magnitude of the error bar. Such a sampling procedure 
is closely related to that used in the finite-temperature 
Lanczos methodi 10 i 14 

As compared with another targeting scheme of finite- 
temperature DMRG)£ the present method has an ad- 
vantage that one does not need to divide the Hilbert 
space with respect to, the z component of total spin, 
but can treat full of the Hilbert space at once. This 
reduces a tedious procedure of numerical simulations sig- 
nificantly. Furthermore, since this technique is of a sim- 
ple extension of the zero temperature DMRG supple- 
mented by the polynomial expansion and random sam- 
pling, momentum-dependent quantities can also be cal- 
culated unlike the case of the transfer-matrix DMRG. 
However, it is not practical to use the present method at 
high temperatures, since the large m is required. We also 
note that, although there is no restriction in principle to 
apply this technique to complicated models with long- 
range interactions, a large number of m is required even 
at low temperatures as is the case of the zero-temperature 
DMRG. 

In summary, a different DMRG technique has been de- 
veloped in order to calculate both static and dynamical 
quantities at low temperatures. This technique is of a 
straightforward extension of a single-target DMRG pro- 
cedure, except that the target state is evaluated by the 
regulated polynomial expansion and a random-sampling 
and averaging procedure are employed for the measure- 
ment of physical quantities. By using the proposed tech- 
nique, static and dynamical quantities in the ID half- 
filled Hubbard chains have been calculated, and it has 
been demonstrated that the technique works well at low 
temperatures. This technique would be useful as one of 
the DMRG techniques at low temperatures. 
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